run_ab<-function(numd,yname,table=1){
pcresults=matrix(NA,ncol=7,nrow=numd)
mrdresults=matrix(NA,ncol=7,nrow=numd)
drdresults=matrix(NA,ncol=7,nrow=numd)	

for (d in 1:numd){
	if (d==1){data <-newdata} else # full sample
	if (d==2){data <-subset(newdata,esse_f>=esse_m&famtype==1&esse_f_nan==0)} else 
	if (d==3){data <-subset(newdata,esse_f<esse_m&famtype==1&esse_f_nan==0)} else
    if (d==4){data <-subset(newdata,famtype==0)} 
      
daty=unlist(data[yname])
datx=data$esse_m
datt=data$qu_m
NT=sum(datx>0)
NC=sum(datx<0)
sss<- (data$Hs/H*N/data$Ns)*(NC/N)
sssT<-(data$Hs/H*N/data$Ns)*(NT/N)
sss[datx>ctff]<-sssT[datx>ctff]
sss2<- (data$Hs/H*N/data$Ns)^2*(NC/N)
sss2T<-(data$Hs/H*N/data$Ns)^2*(NT/N)
sss2[datx>ctff]<-sss2T[datx>ctff]
ygrid=quantile(daty,seq(0.001,0.999,length.out=999))

if (table<=2){
if (d==1){
smooth=4.5
mrdbw.out = rdbwselect(daty,datx, bwselect = "CCT")
bw_mrd=mrdbw.out$bws[1,1]*(NC+NT)^(1/5-1/smooth)	
}
} else
if (table==3){
smooth=4.5
mrdbw.out = rdbwselect(daty,datx, bwselect = "CCT")
bw_mrd=mrdbw.out$bws[1,1]*(NC+NT)^(1/5-1/smooth)	
} else
if (table==4){
smooth=4
if (d==1){
mrdbw.out = rdbwselect(daty,datx, bwselect = "IK")
bw_mrd=mrdbw.out$bws[1,1]*(NC+NT)^(1/5-1/smooth)	
}
}

########## unweighted
if (table==1){
## proportion of complier
results<-MRD_sharp(datt,datx,bw_mrd,Xeval)
pcresults[d,1:3]<-c(bw_mrd,results$estimate,(results$estimate/results$tstat))
## Mean Change in Outcome
results<-MRD_sharp(daty,datx,bw_mrd,Xeval)
mrdresults[d,1:3]<-c(bw_mrd,results$estimate,(results$estimate/results$tstat))
}
########## weighted
## proportion of complier
results<-MRD_sharp(datt,datx,bw_mrd,Xeval,weight=sss)
pcresults[d,4:7]<-c(bw_mrd,results$estimate,(results$estimate/results$tstat),sum(abs(datx)<=bw_mrd))
## Mean Change in Outcome
results<-MRD_sharp(daty,datx,bw_mrd,Xeval,weight=sss)
mrdresults[d,4:7]<-c(bw_mrd,results$estimate,(results$estimate/results$tstat),sum(abs(datx)<=bw_mrd))

# MRD plot
if (table==2){
xlab="Year to Pension Eligibility"
ylab1="Log Equivalized Expenditure"
ylab2="# of Individuals"
main="Mean RD / Histogram"
labels<-data.frame(xlab,ylab1,ylab2,main)
results$bw=bw_mrd
pdf(paste("./APP3/figures/ab_Figure_MRD_CCT",yname,"_d",d,".pdf",sep=""),width=7,height=5)
MRDplot(Xeval,daty,datx,results,labels)
dev.off()
}

## DRD test
if (table>1){
if (table==2|table==3){
smooth=4.5
results<-DRDtest(daty,datx,ygrid,h=bw_mrd,bwsel = "CCT",weight=sss,undersmooth=smooth)
} else
if (table==4){
smooth=4
results<-DRDtest(daty,datx,ygrid,h=bw_mrd,bwsel = "IK",weight=sss,undersmooth=smooth)
}	
drdresults[d,4:7]<-c(results$bw_iter,min(results$ksstat_iter),max(results$ksstat_iter),sum(abs(datx)<=results$bw_iter))
simcv(daty,datx,results$bw_iter,ygrid,weight=sss,weight2=sss2,ccdfC=results$ccdfC_iter,ccdfT=results$ccdfT_iter)
}

## DRD plot
if (table==2){
xlab="Log Equivalized Expenditure"
ylab1="Conditional CDFs"
ylab2="Difference in CCDFs"
main="Distributional RD"	
labels<-data.frame(xlab,ylab1,ylab2,main)
pdf(paste("./APP3/figures/ab_Figure_DRD_CCT",yname,"_d",d,".pdf",sep=""),width=7,height=5)
DRDplot(ygrid,results,labels,iter=1,type=2)
dev.off()
}
}
list(pc=pcresults,mrd=mrdresults,drd=drdresults,pc=pcresults)
}
